////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//---------------------------- Ebrahim Foulaadvand, 11 July 2013 -----------------                        //
//                                                                                                        //
//  The routine "Dalembert" solves the 1D classical wave equation using centred time centered space       //          //
//  CTCS (midpoint leapfrog) scheme. Periodic boundary condition is used. The system length L and the     //
//  wave velocity c are set to one. Initial wave shape is cosine-modulated Gaussian pulse.                //                                                         //
//                                                                                                        //
//                                                                                                        //                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>

using namespace std;

main()
{

double tau=0.004,h,L=1,c=1,a,sigma=0.05,lambda=0.2,k=2*M_PI/lambda,x;

int N=200,i,n,T=1000;


ofstream file0 ("initial profile n=0.plt");  //output file for the profile at timestep n=0.
ofstream file1 ("profile FTCS n=40.plt");    //output file for the profile at timestep n=100.
ofstream file2 ("profile FTCS n=50.plt");    //output file for the profile at timestep n=500.
ofstream file3 ("profile FTCS n=60.plt");    //output file for the profile at timestep n=500.


vector<double> u(N+1,0),unew(N+1,0),uold(N+1,0),f(N+1,0),g(N+1,0);

 // Arrays "u", "unew" and "uold" store the current, updated and previous values
 // of the scaler function u(x,t) at grid points. Array "g" stores the second initial conditions:
 // du(x,0)/dt=g(x). Here we take g(x)=0 for simplicity.

 h=L/N; cout<<"h= "<<h<<endl; cout<<"tau= "<<tau<<endl; //getch();

 a=(c*tau)/(h);

for(i=0;i<=N;i++){

x=-0.5*L+i*h; u[i]=cos(k*x)*exp(-0.5*(x/sigma)*(x/sigma)); g[i]=0.;

file0<<(-0.5*L+i*h)<<"  "<<u[i]<<endl; } //storing the initial conditions data in the data file.



for(i=0;i<=N;i++){ uold[i]=u[i]-g[i]*tau; } // starter condition

unew[0]=2*u[0]-uold[0]+ a*a*(u[1]+u[N]-u[0]); // periodic condition has been used: u[-1] is replaced with u[N]

for(i=1;i<N;i++){ unew[i]=2*u[i]-uold[i] + a*a*(u[i+1]+u[i-1]-2*u[i]); }

unew[N]=2*u[N]-uold[N]+ a*a*(u[0]+u[N-1]-u[N]); // periodic condition has been used: u[N+1] is replaced with u[0]

for(i=0;i<=N;i++){ uold[i]=u[i];  u[i]=unew[i]; }


for (n=1;n<=T;n++){ //cout<<"n= "<<n<<endl; //getch();


// -------------------       CTCS method -------------------------------------


unew[0]=2*u[0]-uold[0]+ a*a*(u[1]+u[N]-2*u[0]); // periodic condition has been used: u[-1] is replaced with u[N]

for(i=1;i<N;i++){ unew[i]=2*u[i]-uold[i] + a*a*(u[i+1]+u[i-1]-2*u[i]); }

unew[N]=2*u[N]-uold[N]+ a*a*(u[0]+u[N-1]-2*u[N]); // periodic condition has been used: u[N+1] is replaced with u[0]


if (n==40)   for(i=0;i<=N;i++){ file1<<(-0.5*L+i*h)<<"  "<<unew[i]<<endl; } //storing in the data file.

if (n==50)   for(i=0;i<=N;i++){ file2<<(-0.5*L+i*h)<<"  "<<unew[i]<<endl; } //storing in the data file.

if (n==60)   for(i=0;i<=N;i++){ file3<<(-0.5*L+i*h)<<"  "<<unew[i]<<endl; } //storing in the data file.


for(i=0;i<=N;i++){ uold[i]=u[i];  u[i]=unew[i];  }

}

cout<<"end." <<endl; getch();

}